Thermodynamic singularities in the entanglement entropy 
at a 2D quantum critical point 
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We study the bipartite entanglement entropy ol the two-dimensional (2D) transverse-field Ising 
model in the thermodynamic limit using series expansion methods. Expansions are developed for 
the Renyi entropy around both the small-field and large-field limits, allowing the separate calcula- 
tion of the entanglement associated with lines and corners at the boundary between sub-systems. 
Series extrapolations are used to extract subleading power laws and logarithmic singularities as the 
quantum critical point is approached. In ID, we find excellent agreement with exact results as well 
as quantum Monte Carlo simulations. In 2D, we find compelling evidence that the entanglement 
at a corner is significantly different from a free boson field theory. These results demonstrate the 
power of the series expansion method for calculating entanglement entropy in interacting systems, a 
fact that will be particularly useful in future searches for exotic quantum criticality in models with 
and without the sign problem. 



Introduction: The study of entanglement properties of 
ground states of one-dimensional (ID) statistical systems 
and free field theories in arbitrary dimensions is a very 
mature field [D-Q. Many exact results have been estab- 
lished, and numerical methods such as the Density Ma- 
trix Renormalization Group (DMRG) enable studies of 
relatively large system sizes in ID [4j. In contrast, the 
study of entanglement properties for ground states of in- 
teracting quantum lattice models in higher dimensions is 
a subject still in its infancy d-Q. In particular, although 
a great potential exists to connect properties of entangle- 
ment to universality at quantum critical points (QCPs) 
[j| , the critical scaling behaviors of very few interacting 
lattice models are known. Ultimately, the study of entan- 
glement entropies may provide unique signatures of novel 
or deconfined QCPs [10j. Yet, much work is required be- 
fore this advance is possible; little is quantitatively known 
about the nature of the singularities and crossovers at a 
QCP as a function of system size and other thermody- 
namic parameters. 

Recent developments in Quantum Monte Carlo meth- 
ods offer a promising avenue for calculating entanglement 
properties of higher dimensional quantum lattice models 
|llL Il2j . Another fruitful approach is the study of en- 
tanglement in suitably parameterized variational wave- 
functions [U, El- DMRG and Matrix Product State 



methods provide other powerful variational approaches 
to study quantum entanglement in higher dimensional 
systems [15l4l7j . However, in contrast to these methods 
that require careful scaling analyses of finite-size lattices, 
series expansions at T = provide a simple yet pow- 
erful alternative approach to studying ground state en- 
tanglement entropy directly in the thermodynamic limit. 
Calculations are carried out order-by-order in perturba- 
tion theory as a power series in some expansion variable 



A, providing a pedagogically transparent introduction to 
the development of entanglement entropy in many-body 
systems. These expansions are typically convergent in- 
side a phase, but become singular as a phase boundary is 
approached. Once the expansions are developed to some 
order (in practice typically of order 10), series extrapo- 
lation methods can be used to approximate the singular 
behavior in entanglement near a QCP. 

Here, we use series expansions to calculate the thermo- 
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where the first sum runs over the nearest-neighbor bonds 
of the square-lattice and the second over its sites. Both 
the limits h = and J = have very simple ground 
states, and series expansions can be separately developed 
in h/J or J/h. At small h one has two ordered ground 
states and the system has spontaneously broken sym- 
metry (called the "ordered" phase). In developing the 
series expansion in h/J, we pick one of the two ground 
states of the system to expand around. The state at 
h = is a simple product state with all the spins pointing 
along the z axis. At large h (the "disordered" phase), one 
also has a simple product ground state, where every spin 
points along the x-axis. Thus both at h = and at J = 
the ground states have no entanglement between any pair 
of sites. A quantum critical point intervenes between the 
ordered and disordered phases, which is known to be in 
the universality class of the 3D classical Ising model [HI . 
Using series expansion, we provide accurate calculation 
of the thermodynamic singularities in the entanglement 
entropy for this universality class, demonstrating in par- 
ticular differences from Gaussian free-field universality. 
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FIG. 1: The infinite square plane is partitioned into four 
quadrants a, b, c and d. The region A could be a half-plane 
such as a U b or a quadrant such as b or c, while the rest of the 
square-lattice forms the region B. For the former partition, 
several low-order clusters that cross the boundary between A 
and B are also shown. 



Series Expansion Methods for Renyi entropies: From 
a computational point of view, the Renyi entropies (l9j 
are particularly convenient measures of bipartite entan- 
glement. If we divide our system into two parts A and 
B, such that each spin belongs to either A or B, then 
the ground state of the full system can be written in the 
local basis as 



(2) 



where a and b refer to basis states for subsystems A and 
B respectively. The Renyi entropies arc defined as: 



1 



1 



ln[Tr(pl)] 



(3) 



where, the trace is over all the states of the subsystem A 
and the reduced density matrix for the subsystem A is 
given by the matrix elements 



{ai\p A \a 2 ) = y~]i>al,bi>a2,b- 



(4) 



In this paper we will focus attention solely on the sec- 
ond Renyi entropy 5*2 . We will divide the infinite system 
into two subsystems such that the subsystem A is either 
a half-plane or a quadrant (See Fig. 1). We begin with 
the case when A is a half-plane. First non-zero terms in 
perturbation theory arise when pairs of spins from across 
the dividing line get entangled. Because the entropy is 
an extensive measure, each such pair contributes equally 
to the sum and it leads to an entropy proportional to the 
length of the boundary. In the next order either a pair 
of spins from one side can be entangled with one spin 
from the other side, or a pair of spins from one side can 
entangle with a pair of spins from the other side. These 
contributions have a natural graphical interpretation in 
terms of clusters that go across the boundary separa ting 
A and B (See Fig. 1). The linked cluster method OH} 
allows one to separate the entanglement that comes from 
a pair of spins versus the additional entanglement that 



comes from a larger cluster of spins. Using the principle 
of inclusion and exclusion one can find the additional en- 
tanglement from a larger cluster (also called the weight 
of the cluster W) by calculating the full entanglement for 
that cluster of spins when the perturbations are turned 
on, and then subtracting from it the weight of all its 
subclusters. 



W(c)=S 2 (c)-J2W(s), 



(5) 



where the sum is over all subclusters of the cluster c. In 
the thermodynamic limit, one can use the translational 
symmetry along the length of the boundary to write the 
entropy per unit length as 



s 2 = S 2 /L = J2w(c d ). 



(6) 



Here the sum is over all translationally distinct clusters 
{cd}- Expanding as a power series in A = h/ J or A = J/h 
one obtains 



S2 = ^p„A" 



(7) 



To obtain the corner terms, we need to consider the 
case where A is a quadrant (See Fig. 1). In fact, by con- 
sidering different choices for the division of the infinite 
lattice into A and B it is possible to completely cancel 
out the line contributions 211. If we calculate the en- 



tanglement entropy for (i) when A is the quadrant (b) 
and (ii) when A is the quadrant (c), then their sum will 
amount to entanglement from two 90 degree corners plus 
two infinite lines that cut across the lattice. The line 
contributions can be subtracted off by subtracting the 
entanglement entropies for the cases, where A is the half 
plane formed by (i) a U b and (ii) a U c. This subtraction 
can be done on a graph by graph basis. Thus series for a 
corner, in the thermodynamic limit, can be expressed in 
terms of graphs that lie at the intersection of two lines. 
This leads to the expansions for the entanglement en- 
tropy at a single corner as 



(8) 



The coefficients p n and q n are calculated to order 14 in 
J/h and order 24 in h/J and provided in supplementary 
material [22]. 

Note that one of the advantages of the series expan- 
sion method is that the line and corner contributions are 
obtained separately. In higher dimensions, entropy asso- 
ciated with each type of manifold, planes, lines, corners 
can also be calculated separately. 

Series analysis: It is clear from the formalism that 
the 'area law' is built into the series expansion method, 
namely that the entanglement entropy scales with the 
boundary 'area' between subsystems. As long as the 
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FIG. 2: Entanglement entropy of the transverse field Ising 
chain, for two edges, obtained by series expansions. For com- 
parison QMC data on finite systems are also shown. In the or- 
dered phase log 2 has been subtracted from the QMC data to 
correspond to the fact that series expansions are done around 
a single ordered ground state. 



perturbation theory converges, the 1 law' continues 
to hold. This is consistent with general arguments that 
gapped phases obey the area law 0] . As one approaches 
a quantum critical point, where the gap goes to zero, the 
series become singular. One can study whether the area 
law continues to hold at the critical point, and the na- 
ture of the critical singularity, by analyzing the limiting 
behavior of the series using extrapolation methods. 

Series extrapolations can deal with convergent or diver- 
gent power-law singularities by using differential approx- 
imants [lH, H3 ■ One represents the function of interest 
/(A) by a differential equation of the form 



df(X) 
dX 



P n (\)f(X) = U j (\). 



(9) 



Here Q m (X), P n (X) and Uj(X) are polynomials of order 
to, n and j, which are obtained such that the differential 
equation correctly reproduces the first m+n+j+2 powers 
in the series expansions for the function /(A). The singu- 
larities of the function arise at A values where Q m {X) = 0. 
If the location of the critical point A c is known by some 
other means, one can put the additional constraint that 
Qm(X c ) = 0. This is called biasing the critical point. One 
can then study the resulting power- law singularity at A c . 
Note that the inhomogeneous term Uj(X) is essential to 
allow a non-zero slowly varying background in addition 
to a power-law singularity. 

In the special case of a log singularity, one can first 
differentiate the function /(A), with respect to A, and 
then represent it by a ratio of polynomials (also called 
Pade approximants), d ^j^ = q" ma ■ The log singularities 
arise for /(A) where Q m {X) = 0. 

First, we discuss the results for the transverse-field 
Ising chain for which closed form expressions for the von 
Neumann entropy and asymptotic expressions for the 
Renyi entropies are available in the literature Our 



goal here is not to study the ID model per sc, but to 
simply sec how well one can estimate the critical proper- 
ties using series expansions of the same length that can 
be done in any dimension. In Fig. 2, the results of sc- 
ries extrapolation are shown. As a comparison, we also 
show finite-size data from Quantum Monte Carlo (QMC) 
on an L length chain with periodic boundary conditions, 
where A and B are both of length L/2. The QMC was 
performed using a T = projector method with cluster 
updates |25[ , adapted to calculate Renyi entropy via the 
Swap operator 11 1 on a replicated system . The com- 
parison shows that both the series expansion and QMC 
results are very accurate, at least until one gets very close 
to the critical point, where finite size effects become large 
and the QMC data drops away from the series extrapo- 
lation curve. 

In ID, the boundary 'area' between subsystems is just 
a point or a corner. At the critical point this corner con- 
tribution diverges logarithmically leading to a breakdown 
of the 'area law'. The second Renyi entropy associated 
with a single boundary is given asymptotically close to 
the critical point by the expression 0,1^1, 



So 



s 



log£, 



(10) 



where the central charge c = | for this model and the 
correlation length £ diverges as 1/|1 — A| as A approaches 
unity. This means that the coefficient of the logarith- 
mic singularity in log |1 — A| should equal —0.0625. Se- 
ries extrapolations (with the length of series compara- 
ble to what is calculated in 2D) give different answers 
from the expansions in h/J and J/h. From one side 
one estimates the coefficient of the log singularity to 
be —0.053(1) where as from the other side one obtains 
—0.077(2). The internal consistency of the approximants 
leads to an unusually small estimate of the systematic 
uncertainty. In fact, one finds that with increasing order 
both terms are changing in the right direction but only 
by about 0.001 and 0.002 respectively in each order. If 
we make the reasonable assumption that the coefficient 
must be the same from both sides and thus average the 
two answers, one would obtain —0.065(12), which gives 
a stringent limit on the uncertainty in the calculations. 

We now turn to the 2D transverse-field Ising model. 
The line and corner terms for the entanglement entropy 
of the 2D transverse- field Ising model are shown in Fig. 3. 
In this case, the critical point is biased to the value of 
h/J 



27 



3.044, a value determined previously [1 
Clearly the line term has a sharp peak at the critical 
point, but it does not diverge, implying that the area law 
continues to hold at the critical point. One expects the 
power-law singularity to be of the form | A — A c | " , where v 
is the critical exponent for the divergence of correlation 
length in the 3D classical Ising model. Our series extrap- 
olations lead to estimates for v = 0.60(2) and entropy per 
unit length at the critical point of s\ = 0.0324(3) from 
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FIG. 3: 'Area law' term S2 and corner term C2 of the en- 
tanglement entropy of the 2D transverse-field Ising model ob- 
tained by series expansions in the variables h/J and J/h. 
In each case two approximants with critical point biased at 
h/J = 3.044 are shown. 



one side and v — 0.66(3) and s 2 = 0.0350(3) from the 
other side. Averaging these we get, v = 0.63(3) and 
s 2 = 0.337(13). These values are clearly consistent with 
the known value of v = 0.629(2) [H] and recent QMC 
estimate of s\ = 0.0332(4) 

The corner term is biased to have a log singularity. We 
estimate the coefficient of the logarithm to be 0.0059(3) 
from one side and 0.0077(1) from the other side. Aver- 
aging the two, we get 0.0068(9). One can convert the 
logarithm in the variable |A C — A into log(£) by dividing 
by — v 1 where v is the correlation length exponent. Thus, 
we estimate, asymptotically, for a single corner 



c 2 = (-0.011 ±0.001) log(£). 



(11) 



These results are clearly distinct from the free field theory 
result of Casini and Huerta who obtain c 2 = —0.0062 
[l^. The QMC in Ref. quotes a value of 4c 2 = 
—0.03 ± 0.01, with large uncertainties that can not be 
distinguished from free field theory. 

Discussions: We have shown that series expansions 
can be used to calculate thermodynamic singularities in 
the entanglement entropy at 2D quantum critical points 
(QCPs) with about 10 percent accuracy. We have pro- 
vided compelling evidence that the entanglement entropy 
produced at a corner in the boundary between subregions 
in the 2D transverse field Ising model QCP is different 
from that predicted in a free boson field theory [3(| . In- 
deed, the transverse field Ising model QCP is the classical 
3D Ising universality class, which is distinct from the free 
(Gaussian) universality class. Currently there are no the- 
ory predictions for the corner log for the transverse field 
Ising model. 

Through our calculations of c 2 (and confirmed by our 
accurate estimate of the exponent v), we have demon- 
strated that series expansions already suffice to distin- 
guish between different universality classes [3(1 31 1. It 



lattices to further consolidate the notion of universality. 
Given that the numerical values of critical parameters 
are largely unknown, comparison between series expan- 
sions and QMC data, where available, would be most 
useful. Series expansions can also be developed in higher 
than two dimensions and also for other Renyi indices n. 
These can help address questions related to upper critical 
dimensionality, boundary correlation functions and pos- 
sible singularities as a function of the Renyi index n (3). 
For example, how does the logarithmic singularity associ- 
ated with the corners depend on dimensionality, and docs 
it have a simple limit above the upper critical dimension? 

A weakness of the series method, as presented here, 
is the inability to study topological entanglement en- 
tropy. The calculations discussed here are all related to 
the boundary between subsystems and in finite orders 
of perturbation theory only the degrees of freedom at fi- 
nite distance from the boundary get entangled. This can 
not address topological entanglement, which is inherently 
long-ranged. It would be interesting to explore the pos- 
sibility of addressing this through an approach involving 
degenerate perturbation theory (32| . 

From a computational point of view, series expansion 
methods may be particularly useful in studying interact- 
ing Fermion models and frustrated spin models. Quan- 
tum critical points are, in principle, accessible to high 
temperature series expansions, which might provide a 
useful route to studying t — J and Hubbard models. At 
T = 0, one should be able to look for exotic critical points 
at the boundary of magnetically ordered phases, or at the 
boundary between ordered phases and gap ped spin liq- 
uids [13, HH- Indeed, several recent works |l5l. [la] have 
argued that a spin liquid state may arise in the frustrated 
J\ — J 2 square-lattice Heisenberg model. Investigations 
of the entanglement scaling at critical points contained 
in this model will be pursued in future. 
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would be useful to study a range of models on different 
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Supplementary Material for "Thermodynamic 
singularities in the entanglement entropy at a 2D 
quantum critical point" by R. R. P. Singh, R. G. 
Melko and J. Oitmaa 

This supplement provides the reader with further tech- 
nical data in support of the main part of the paper. First 
we present tables of series coefficients and then tables 
showing estimates of the log singularity 



n 


Pn 


2 


0.125 


4 


0.4921875 


6 


2.38053385 


8 


14.0315145 


10 


90.7878208 


12 


625.42033 


14 


4501.17322 



TABLE I: Series expansion coefficients for entanglement en- 
tropy per unit length S2 (p n in Eq. 7 in the paper) for 2D 
Transverse Field Ising model in the variable J/h. 



n 


Pn 


4 


5.425349E-05 


6 


2.9669904E-06 


8 


3.61269111E-07 


10 


2.23023576E-08 


12 


2.64005223E-09 


14 


1.64460642E-10 


16 


2.27121391E-11 


18 


1.43038774E-12 


20 


1.84792391E-13 


22 


1.47160415E-14 


24 


1.56354826E-15 



TABLE II: Series expansion coefficients for entanglement en- 
tropy per unit length S2 (p n in Eq. 7 in the paper) for 2D 
Transverse Field Ising model in the variable h/J. 



n 


q n 


4 


-0.21875 


6 


-1.578125 


8 


-12.1454976 


10 


-94.8888162 


12 


-757.353415 


14 


-6151.86971 



TABLE III: Series expansion coefficients for entanglement en- 
tropy of a corner C2 (q n in Eq. 8 in the paper) for 2D Trans- 
verse Field Ising model in the variable J/h. 



6 



n 


q-n 


6 


-4 


23855252E- 


07 


8 


-7 


99438098E- 


08 


10 


-9 


17332192E- 


09 


12 


-8 


55136841E- 


10 


14 


-1 


02753121E- 


10 


16 


-9 


44954642E- 


12 


18 


-1 


09399979E 


12 


20 


-9 


94186394E- 


14 


22 


-1 


13706509E 


14 


24 


-1 


08159628E 


15 



TABLE IV: Series expansion coefficients for entanglement en- 
tropy of a corner a (q n in Eq. 8 in the paper) for 2D Trans- 
verse Field Ising model in the variable h/J. 



{J/h)l 


amplitude 


m 


n 


0.107922179 


0.00762850877 


2 


4 


0.107922179 


0.00772136388 


2 


5 


0.107922179 


0.00754811631 


3 


3 


0.107922179 


0.00775324771 


3 


4 


0.107922179 


0.00761989088 


4 


2 


0.107922179 


0.00786590457 


4 


3 


0.107922179 


0.00777476393 


5 


2 



TABLE V: Pade approximant estimates for amplitude of log- 
divergence for 2D Transverse Field Ising model in the variable 
(J/h) 2 . The first column is the biased location of the critical 
point. The second is the amplitude for the log divergence. 
The third and fourth give m n values for the Pade approxi- 
mants. 



{h/J)l 


amplitude 


m 


n 


9 


265936 


0.00619029914 


3 


6 


9 


265936 


0.00590459618 


3 


7 


9 


265936 


0.00573168376 


3 


8 


9 


265936 


0.00698086931 


4 


5 


9 


265936 


0.00586117615 


4 


6 


9 


265936 


0.0053144883 


4 


7 


9 


265936 


0.00588231756 


4 


8 


9 


265936 


0.00578324783 


5 


4 


9 


265936 


0.00624218823 


5 


5 


9 


265936 


0.00560957049 


5 


6 


9 


265936 


0.00584200216 


5 


7 


9 


265936 


0.00608591271 


6 


4 


9 


265936 


0.00604662478 


() 


5 


9 


265936 


0.00613351483 


6 


6 


9 


265936 


0.00604144926 


7 


4 


9 


265936 


0.00607802345 


7 


5 


9 


265936 


0.00624832351 


8 


4 



TABLE VI: Pade approximant estimates for amplitude of log- 
divergence for 2D Transverse Field Ising model in the variable 
(h/J) 2 . The first column is the biased location of the critical 
point. The second is the amplitude for the log divergence. 
The third and fourth give m n values for the Pade approxi- 
mants. 



